library(tidyverse)
library(ComplexHeatmap)
library(circlize)
library(pheatmap)
library(corrplot)
library(ggplot2)
library(reshape2)
library(grid)
library(dplyr)


# Read the tables before running the R script.

#------------------- start -----------------------



#--------------------------------------------------
#----------------- Xenium level1 -----------------
#--------------------------------------------------

i = 1
df <- t(df_filtered_list[[i]])

xenium_order = c(
  "Xe_anterior forebrain",
  "Xe_infundibular organ (IO)",
  "Xe_dorsal layer",
  "Xe_cells of Joseph",
  "Xe_dorsal ventricular cavity",
  "Xe_translumenal cells I (anterior)",
  "Xe_boundary cells I",
  "Xe_boundary cells II",
  "Xe_translumenal cells II (posterior-upper)",
  "Xe_translumenal cells III (posterior-lower)",
  "Xe_central canal",
  "Xe_nucleus of Rohde (nRo)",
  "Xe_migrated cells",
  "Xe_others"
)

snRNA_order = c(
  "Bf_21", "Bf_31", "Bf_33", "Bf_16", "Bf_30", "Bf_35", "Bf_11", "Bf_38", "Bf_0", "Bf_23",
  "Bf_6", "Bf_17", "Bf_28", "Bf_15", "Bf_18", "Bf_34",
  "Bf_13", "Bf_3", "Bf_39", "Bf_9", "Bf_36", "Bf_22", "Bf_7", "Bf_1",
  "Bf_2", "Bf_8", "Bf_19", "Bf_26", "Bf_4", "Bf_12", "Bf_24", "Bf_25", "Bf_5",
  "Bf_37", "Bf_20", "Bf_27", "Bf_10", "Bf_14", "Bf_32", "Bf_29"
)

level1_top_anno = HeatmapAnnotation(
  type = factor(colnames(df), levels = xenium_order),
  col = list(type = c(
    "Xe_anterior forebrain" = "#CF8EB1",
    "Xe_boundary cells I" = "#CF8E4D",
    "Xe_boundary cells II" = "#A65628",                      
    "Xe_central canal" = "#549F6C",
    "Xe_dorsal layer" = "#FF7F00",
    "Xe_dorsal ventricular cavity" = "#F1E330",
    "Xe_infundibular organ (IO)" = "#E41A1C",
    "Xe_cells of Joseph" = "#C1862A",
    "Xe_translumenal cells I (anterior)" = "#FECB20",
    "Xe_translumenal cells II (posterior-upper)" = "#984EA3",
    "Xe_translumenal cells III (posterior-lower)" = "#728A7A",
    "Xe_migrated cells" = "#BFD1C4",
    "Xe_nucleus of Rohde (nRo)" = "#377EB8",
    "Xe_others" = "#9C3EFF"
  )),
  annotation_name_gp = gpar(fontsize = 8),
  simple_anno_size = unit(2.5, "mm"),
  show_annotation_name = FALSE,
  show_legend = T
)
level1_row_anno = rowAnnotation(
  type = c( #Bf_0, Bf_1, Bf_10, Bf_11, Bf_12...
    "Non-neuron", "Neuron", "Neuron", "Non-neuron", "Non-neuron", "Neuron", "Neuron", "Neuron",
    "Neuron", "Neuron", "Neuron", "Non-neuron", "Non-neuron", "Neuron", "Neuron", "Neuron",
    "Non-neuron", "Neuron", "Neuron", "Neuron", "Neuron", "Neuron", "Neuron", "Neuron",
    "Neuron", "Neuron", "Neuron", "Neuron", "Neuron", "Neuron", "Neuron", "Neuron",
    "Non-neuron", "Neuron", "Non-neuron", "Neuron", "Neuron", "Neuron", "Neuron", "Neuron"
  ),
  col = list(type = c(  
    "Neuron" = "#dab1da",
    "Non-neuron" = "lightgreen"
  )),
  annotation_name_gp = gpar(fontsize = 8),
  simple_anno_size = unit(2.5, "mm"),
  show_annotation_name = FALSE,
  show_legend = T
)



# generate heatmap
ht_list = Heatmap(df,
        top_annotation = level1_top_anno,
        right_annotation = level1_row_anno,
        col = colorRamp2(seq(0, 1, length.out = 25), c("#fff",rev(hcl.colors(24, "Reds 2")))),
        name = name[[i]],
        column_names_rot = 45,
        column_order = xenium_order,
        column_names_side = "top",
        row_order = snRNA_order,
        row_names_side = "right",
        show_row_dend = TRUE,
        width = ncol(df) * unit(8, "mm"),
        height = nrow(df) * unit(4.5, "mm"),
        rect_gp = gpar(col = "#aaa", lwd = 1),
        row_names_gp = gpar(fontsize = 12),
        column_names_gp = gpar(fontsize = 12),
        show_heatmap_legend = TRUE,
        border = TRUE,
        cell_fun = function(j, i, x, y, width, height, fill) {
          if (df[i, j] > 0.2) {
            grid.points(x, y, pch = 16, size = unit(1.5, "mm"), gp = gpar(col = "black"))
          }
        }
)


ht_drawn = draw(ht_list)


w = ComplexHeatmap:::width(ht_drawn)
h = ComplexHeatmap:::height(ht_drawn)

w_inch = convertWidth(w, "inch", valueOnly = TRUE)
h_inch = convertHeight(h, "inch", valueOnly = TRUE)

pdf("./Bf snRNAseq_r1 vs Xenium_Level1.pdf", width = w_inch, height = h_inch)
draw(ht_list)
dev.off()




#--------------------------------------------------
#----------------- Xenium cluster -----------------
#--------------------------------------------------

i = 2
df <- t(df_filtered_list[[i]])

Xe_cluster_order = c(
  "Xe_11", "Xe_3", "Xe_33", "Xe_16", "Xe_29", "Xe_26", "Xe_32", "Xe_2",
  "Xe_17", "Xe_25", "Xe_24", "Xe_6", "Xe_18", "Xe_8", "Xe_9", "Xe_21",
  "Xe_13", "Xe_19", "Xe_5", "Xe_4", "Xe_35", "Xe_1", "Xe_0", "Xe_10",
  "Xe_22", "Xe_7", "Xe_15", "Xe_20", "Xe_14", "Xe_12", "Xe_27", "Xe_28",
  "Xe_23", "Xe_31", "Xe_30", "Xe_34"
)

top = c(
  "Xe_translumenal cells III (posterior-lower)",
  "Xe_translumenal cells III (posterior-lower)",
  "Xe_central canal",
  "Xe_anterior forebrain",
  "Xe_nucleus of Rohde (nRo)",
  "Xe_boundary cells I",
  "Xe_nucleus of Rohde (nRo)",
  "Xe_central canal",
  "Xe_anterior forebrain",
  "Xe_cells of Joseph",
  "Xe_dorsal ventricular cavity",
  "Xe_boundary cells II",
  "Xe_dorsal layer",
  "Xe_central canal",
  "Xe_translumenal cells I (anterior)",
  "Xe_central canal",
  "Xe_migrated cells",
  "Xe_cells of Joseph",
  "Xe_cells of Joseph",
  "Xe_anterior forebrain",
  "Xe_migrated cells",
  "Xe_migrated cells",
  "Xe_anterior forebrain",
  "Xe_anterior forebrain",
  "Xe_others",
  "Xe_migrated cells",
  "Xe_infundibular organ (IO)",
  "Xe_anterior forebrain",
  "Xe_others",
  "Xe_translumenal cells II (posterior-upper)",
  "Xe_translumenal cells II (posterior-upper)",
  "Xe_translumenal cells II (posterior-upper)",
  "Xe_dorsal ventricular cavity",
  "Xe_central canal",
  "Xe_translumenal cells I (anterior)",
  "Xe_translumenal cells I (anterior)"
)

top_order = c(
  "Xe_anterior forebrain",
  "Xe_infundibular organ (IO)",
  "Xe_dorsal layer",
  "Xe_cells of Joseph",
  "Xe_dorsal ventricular cavity",
  "Xe_translumenal cells I (anterior)",
  "Xe_boundary cells I",
  "Xe_boundary cells II",
  "Xe_translumenal cells II (posterior-upper)",
  "Xe_translumenal cells III (posterior-lower)",
  "Xe_central canal",
  "Xe_nucleus of Rohde (nRo)",
  "Xe_migrated cells",
  "Xe_others"
)


top_anno = HeatmapAnnotation(
  type = factor(top, levels = top_order),
  col = list(type = c(
    "Xe_anterior forebrain" = "#CF8EB1",
    "Xe_boundary cells I" = "#CF8E4D",
    "Xe_boundary cells II" = "#A65628",                      
    "Xe_central canal" = "#549F6C",
    "Xe_dorsal layer" = "#FF7F00",
    "Xe_dorsal ventricular cavity" = "#F1E330",
    "Xe_infundibular organ (IO)" = "#E41A1C",
    "Xe_cells of Joseph" = "#C1862A",
    "Xe_translumenal cells I (anterior)" = "#FECB20",
    "Xe_translumenal cells II (posterior-upper)" = "#984EA3",
    "Xe_translumenal cells III (posterior-lower)" = "#728A7A",
    "Xe_migrated cells" = "#BFD1C4",
    "Xe_nucleus of Rohde (nRo)" = "#377EB8",
    "Xe_others" = "#9C3EFF"
  )
  ),
  annotation_name_gp = gpar(fontsize = 8),
  simple_anno_size = unit(2.5, "mm"),
  show_annotation_name = FALSE,
  show_legend = T
)


# generate heatmap
ht_list = Heatmap(df,
        top_annotation = top_anno,
        right_annotation = level1_row_anno,
        col = colorRamp2(seq(0, 1, length.out = 25), c("#fff",rev(hcl.colors(24, "Reds 2")))),
        column_names_rot = 45,
        column_order = Xe_cluster_order,
        column_names_side = "top",
        row_names_side = "right",
        row_order = snRNA_order,
        show_row_dend = TRUE,
        width = ncol(df) * unit(6, "mm"),
        height = nrow(df) * unit(4.5, "mm"),
        rect_gp = gpar(col = "#aaa", lwd = 1),
        row_names_gp = gpar(fontsize = 12),
        column_names_gp = gpar(fontsize = 12),
        show_heatmap_legend = TRUE,
        cell_fun = function(j, i, x, y, width, height, fill) {
          if (df[i, j] > 0.2) {
            grid.points(x, y, pch = 16, size = unit(1.5, "mm"), gp = gpar(col = "black"))
          }
        },
        border = TRUE)


ht_drawn = draw(ht_list)


w = ComplexHeatmap:::width(ht_drawn)
h = ComplexHeatmap:::height(ht_drawn)

w_inch = convertWidth(w, "inch", valueOnly = TRUE)
h_inch = convertHeight(h, "inch", valueOnly = TRUE)

pdf("./Bf snRNAseq_r1 vs Xenium_cluster.pdf", width = w_inch, height = h_inch)
draw(ht_list)
dev.off()
